% Measure intervals between strobed word digital outputs.
% @param dOutObject an instance of a dotsAllDOutObjects subclass
% @param getExternalFunction a function that will return the
% interval between the last two strobed words as "seen" from an external
% data collection system (such as Plexon, a USB device, etc.).  The
% external data collection system must be connected to the outputs used by
% @a dOutObject.
% @param waitFunction optional function that will wait or sleep for a given
% duration.  Defaults to Matlab's builtin pause().
% @param clockFunction optional function that will return the current time
% as a number.  Defaults to the clockFunction shared value in.
% dotsTheSwitchboard.
% @param words optional words to ouput from @a dOutObject.  Defaults to
% [0:255].
% @param intervals optional delays to wait between repetition of each word
% in @a words.  Defaults to linspace(0, 1, 10).
% @details
% <b>The Test</b>
% @details
% Uses @a dOutObject to output m x n pairs of strobed words, separated
% by various delays, where m is the number of delays and n is the number of
% strobed words.  The delay values are specified in @a intervals and the
% delays carried out by @a waitFunction.  The strobed word values are
% specified in @a words.  Each word value is output twice per pair.
% @details
% Outputs each pair of strobed words with the following sequence:
%   - uses @a dOutObject to output the first of a word pair and gets a
%   "before timestamp"
%   - uses @a clockFunction to get a "before clock time"
%   - uses @a waitFunction to delay for a given "nominal interval"
%   - uses @a dOutObject to output the last of a word pair and gets an
%   "after timestamp"
%   - uses @a clockFunction to get an "after clock time"
%   - evaluates @a getExternalFunction to access the interval between
%   strobed words as "seen" from an external data collection system.
%   .
% Thus, each "nominal interval" can be compared to three other measurements
% of the same interval: the interval between @a dOutObject timestamps, the
% interval between @a clockFunction clock times, and the interval between
% strobed words as seen form an external data collections system.
% @details
% <b>The Results</b>
% @details
% In the ideal, all three interval measurements would match the nominal
% interval.  But in practice, each function call takes a finite, random
% time to complete.  Furthermore, the a local clock and the clock of an
% external data collection system mgith count at different rates.  So the
% measured intervals may differ from nominal, and from one another.
% @details
% benchmarkDOutStrobedWords makes three plots to compare interval
% measurements.  The top plot assesses basic agreement among the nominal
% intervals and the three measured intervals.  It's internded as a raw
% "smell test" of the raw measured intervals: its time scale is coarse and
% each data point may hide a cluster of data points.
% @details
% The bottom two plots focus on how well the intervals between @a
% dOutObject timestamps agree with the intervals as as "seen" from the
% external data collection system.
% @details
% For the bottom left plot, intervals from @a dOutObject and @a
% getExternalFunction are grouped by nominal interval.  The mean of each
% group is subtracted, so that all the intervals cluster near 0.  The
% mean-subtracted @a dOutObject and @a getExternalFunction intervals are
% plotted against each other.  If the measurements from each source, agree,
% all the points will fall near the unity line (i.e. x=y).
% @details
% For the bottom right plot, intervals from @a dOutObject and @a
% getExternalFunction are each transformed by simple linear regression into
% the time frame of @a waitFunction (i.e. the nominal delay).  See
% estimateClockDriftAndOffset() for more about transforming time values as
% "seen" from different clocks).  This transformation accounts the mutual
% offset and drift that is to be expected between clocks, in this case
% between the the clock used by @a waitFunction and each clock used by @a
% dOutObject or the clock from the external data collection system.
% @details
% The residual "jitteriness", not not accounterd for by the linear
% transaformation, is plotted for the @a dOutObject and @a
% getExternalFunction intervals.  Again, the points should fall along the
% unity line.
% @details
% Returns a struct with 5 fields:
%   - @b words the given @a words, or defaults
%   - @b intervals the given delay @a intervals, or defaults
%   - @b dOutResiduals m x n matrix of jitter residuals, as "seen" from the
%   timestamps returned from @a dOutObject.  m is the number of @a
%   intervals and n is the number of @a words
%   - @b getExternalResiguals mxn matrix of jitter residuals, as "seen"
%   from the timestamps returned by getExternalFunction
%   - @b clockResiduals mxn matrix of jitter residuals, as "seen"
%   from the given @clockFunction
%   .
% @details
% benchmarkDOutStrobedWords leaves some details of how to run the test up
% to its caller.  In particular, it does not "know" what type of hardware
% dOutObject uses or how @a dOutObject implements the dotsAllDOutObjects
% interface.  It also does not "know" to what external data collection
% system @a dOutObject is connected or what hardware or software are
% involved.  It only expectes that when called @a getExternalFunction will
% return the time interval between the two most recent strobed word
% outputs.
% %details
% See benchmarkSnowDots for a usage example of benchmarkDOutStrobedWords.
%
% @ingroup utilities
function data = benchmarkDOutStrobedWords(dOutObject, ...
    getExternalFunction, waitFunction, clockFunction, words, intervals)

data = [];

% validate the object that will send digital outputs
if nargin < 1 || ...
        ~isa(dOutObject, 'dotsAllDOutObjects') ...
        || ~dOutObject.isAvailable
    disp(sprintf('%s needs a working dotsAllDOutObjects object', mfilename))
    return
end

% use a Matlab-builtin default wait function
if nargin < 3 || isempty(waitFunction)
    waitFunction = @pause;
end

% use a platform-correct default wait function
%   which should be set with dotsTheMachineConfiguration
if nargin < 4 || isempty(clockFunction)
    clockFunction = ...
        dotsTheSwitchboard.getSharedPropertyValue('clockFunction');
end

% use 8-bit words by default
if nargin < 5 || isempty(words)
    words = 0:255;
end

% use 10 different intervals by default
if nargin < 6 || isempty(intervals)
    intervals = linspace(0, 1, 10);
end

disp(sprintf('%s: expected test duration: %.2f by %s', ...
    mfilename, sum(intervals)*numel(words), func2str(waitFunction)));

n = numel(words);
m = numel(intervals);
dOutIntervals = zeros(m,n);
clockIntervals = zeros(m,n);
getExternalIntervals = zeros(m,n);

% send word pairs and collect timestamps
%   duplicate the first word to make sure functions are loaded
for jj = [1 1:n]
    for ii = [1 1:m]
        dOutPre = dOutObject.sendStrobedWord(words(jj), 0);
        clockPre = feval(clockFunction);
        feval(waitFunction, intervals(ii));
        clockPost = feval(clockFunction);
        dOutPost = dOutObject.sendStrobedWord(words(jj), 0);
        
        dOutIntervals(ii,jj) = dOutPost - dOutPre;
        clockIntervals(ii,jj) = clockPost - clockPre;
        getExternalIntervals(ii,jj) = feval(getExternalFunction);
    end
end

% organize data to be returned
data.words = words;
data.intervals = intervals;
data.waitFunction = waitFunction;
data.clockFunction = clockFunction;

nominalIntervals = repmat(intervals', 1, n);

% For each interval measurement, compute the simple regression to map
% measured intervals into the frame of reference of the nominal intervals
% and compute the residuals jitter not accounted for by the regression
dOutToNominal = estimateClockDriftAndOffset( ...
    dOutIntervals, nominalIntervals);
dOutFit = applyClockDriftAndOffset( ...
    dOutIntervals, dOutToNominal);
data.dOutResiduals = nominalIntervals - dOutFit;
data.dOutIntervals = dOutIntervals;

clockToNominal = estimateClockDriftAndOffset( ...
    clockIntervals, nominalIntervals);
clockFit = applyClockDriftAndOffset( ...
    clockIntervals, clockToNominal);
data.clockResiduals = nominalIntervals - clockFit;
data.clockIntervals = clockIntervals;

getExternalToNominal = estimateClockDriftAndOffset( ...
    getExternalIntervals, nominalIntervals);
getExternalFit = applyClockDriftAndOffset( ...
    getExternalIntervals, getExternalToNominal);
data.getExternalResiduals = nominalIntervals - getExternalFit;
data.getExternalIntervals = getExternalIntervals;

%% plot strobed word data
f = figure(100);
clf(f);
set(f, 'Name', mfilename, 'NumberTitle', 'off');

% the top axes show raw data
spacer = min(diff(intervals));
rawAx = subplot(2,1,1, ...
    'Parent', f, ...
    'XLim', [min(intervals)-spacer, max(intervals)+spacer], ...
    'XTick', unique(sort(intervals)), ...
    'XGrid', 'on', ...
    'YLim', [min(intervals)-spacer, max(intervals)+spacer], ...
    'YTick', unique(sort(intervals)), ...
    'YGrid', 'on');
title(rawAx, 'strobe pair delay intervals')
ylabel(rawAx, 'measured interval')
xlabel(rawAx, ...
    sprintf('nominal interval from %s', func2str(waitFunction)))

% the bottom left axes show mean-subtracted data for local and external
% estimates of the interval between digital outputs
lims = [-.002 .002];
ticks = -.002:.001:.002;
subtractAx = subplot(2,2,3, ...
    'XLim', lims, ...
    'XTick', ticks, ...
    'XGrid', 'on', ...
    'YLim', lims, ...
    'YTick', ticks, ...
    'YGrid', 'on', ...
    'Parent', f);
title(subtractAx, 'mean-subtracted')
xlabel(subtractAx, 'dOutIntervals')
ylabel(subtractAx, 'getExternalIntervals')

% the bottom left axes show regression residuals for local and external
% estimates of the interval between digital outputs
residualAx = subplot(2,2,4, ...
    'XLim', lims, ...
    'XTick', ticks, ...
    'XGrid', 'on', ...
    'YLim', lims, ...
    'YTick', ticks, ...
    'YGrid', 'on', ...
    'Parent', f);
title(residualAx, 'fit-to-nominal residual')
xlabel(residualAx, 'dOutIntervals')
ylabel(residualAx, 'getExternalIntervals')

intervalNames = { ...
    'clockIntervals', 'dOutIntervals', 'getExternalIntervals'};

nNames = numel(intervalNames);
cols = lines(nNames);
offOn = {'off', 'on'};
for ii = 1:m
    for jj = 1:nNames
        intervalTimes = data.(intervalNames{jj});
        placement = zeros(1,n)+(intervals(ii)+((jj-nNames/2)*(spacer/10)));
        line(placement, intervalTimes(ii,:), ...
            'Parent', rawAx, ...
            'HandleVisibility', offOn{1+(ii==1)}, ...
            'Marker', '+', ...
            'LineStyle', 'none', ...
            'Color', cols(jj,:));
    end
    
    line(data.dOutIntervals(ii,:) - mean(data.dOutIntervals(ii,:)), ...
        data.getExternalIntervals(ii,:) - mean(data.getExternalIntervals(ii,:)), ...
        'Parent', subtractAx, ...
        'Marker', '+', ...
        'LineStyle', 'none', ...
        'Color', [0 0 0]);
    
    line(data.dOutResiduals(ii,:), data.getExternalResiduals(ii,:), ...
        'Parent', residualAx, ...
        'Marker', '+', ...
        'LineStyle', 'none', ...
        'Color', [0 0 0]);
end
legend(rawAx, intervalNames{:}, 'Location', 'best')